Freely decaying weak turbulence for sea surface gravity waves 
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We study numerically the generation of power laws in the framework of weak turbulence theory 
for surface gravity waves in deep water. Starting from a random wave field, we let the system evolve 
numerically according to the nonlinear Euler equations for gravity waves in infinitely deep water. In 
agreement with the theory of Zakharov and Filonenko, we find the formation of a power spectrum 
characterized by a power law of the form of |k| -2 ' 5 . 



After the pioneering work by Kolmogorov ffl on the 
equilibrium range in the spectrum of an homogeneous 
and isotropic turbulent flow, there have been a number 
of studies on cascade processes in many other fields of 
classical physics such as plasma physics, magnetohydro- 
dynamics and ocean waves. For surface gravity waves the 
first seminal theoretical work was done by O.M. Phillips 
in 1958, Using dimensional arguments, he argued 
that the frequency spectrum in the inertial range was of 
the form F(u>) = ag 2 uj~ 5 , where a was supposed to be 
an absolute constant and g is gravity. Even though in 
the introduction of Phillips' paper it was stated that "a 
necessary condition for the equilibrium range over a cer- 
tain part of the spectrum is the appreciable non-linear 
interactions among these wave- numbers" (from g]), his 
arguments were based on the geometrical features of the 
free surface elevation. One of his basic assumptions was 
that the only variable of interest was gravity, while the 
friction velocity, it*, was not supposed to be involved in 
the spectral relation, limiting the possibility for a correct 
dimensional analysis. 

Some years later Zakharov and Filonenko j|] estab- 
lished that in infinite water the direct cascade should 
produce a power spectrum of the surface elevation of the 
form P(|k|) ~ |k|~ 2 - 5 that corresponds, using the linear 
dispersion relation in infinite depth, to an w -4 frequency 
power spectrum: the result was found as an exact solu- 
tion of the kinetic wave equation (see Q ) . The theory de- 
veloped is known as "weak" or "wave turbulence" and has 
many important applications in different fields of physics 
such as hydrodynamics, plasma physics, nonlinear optics, 
solid state physics, etc. see ||. It is called weak turbu- 
lence because it deals with resonant interactions among 



small-amplitude waves. Thus, contrary to fully devel- 
oped turbulence, it leads to explicit analytical solutions 
provided some assumptions are made. The first exper- 
imental support of the theory for surface gravity waves 
was made by Toba Q who was completely unaware of 
the paper by Zakharov and Filonenko. He reformulated 
the Phillips' equilibrium range law in the following way: 
F(uj) — (3gu^uj~ A , where (3 should now be a universal di- 
mensionless constant. After the work by Toba, successive 
experimental observation of the u>~ 4 law have been made 
by a number of authors, see for example [f7|dTo|. 

Even though there is a consensus on this result, it must 
be stressed that so far the verification of the theory has 
never been established from first principles and more- 
over the mechanisms that lead to the power law w~ 4 
are not universally recognized: geometrical aspects re- 
lated to wave breaking, without invoking the nonlinear 
wave-wave interaction mechanism, are still retained by 
many oceanographers as fundamental for generating an 

power law. Confirmation of the Zakharov-Filoncnko 
solution to the kinetic equation has being given through 
numerical simulations of the kinetic wave equation itself 

|Q, solving exactly the so called S n i term. Never- 
theless, it must be underlined that the kinetic equation is 
derived from the primitive equations of motion under a 
number of hypotheses (see for example fl5)| ), therefore it 
cannot be concluded a priori that power law solutions of 
the kinetic equation are also shared by the fully nonlinear 
wave equations. 

One way to verify weak turbulence theory is to per- 
form direct numerical simulations of the primitive equa- 
tions of motion. The numerical confirmation of the the- 
ory for gravity waves propagating on a surface has not 
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been an easy task (for capillary waves see |1J], for one 
dimensional wave turbulence see [tL5Ul6fl ), basically be- 
cause of the intrinsic difficulties of the computation of 
the boundary conditions. Different numerical approaches 
have been used for integrating the fully nonlinear sur- 
face gravity waves equations (see (T^J for a review) . The 
numerical methods based on volume formulations show 
very interesting results, in particular they are capable 
of modeling in a quite appropriate way wave breaking. 
Unfortunately they have the disadvantage that they re- 
quire large computational resources, and therefore are 
not suitable for long time numerical simulations. For ir- 
rotational and inviscid flows boundary formulations are 
usually preferred: only the surface is discretized reduc- 
ing the dimension of computation (from three to two). 
The Higher-Order Spectral Methods (HOS), indeed the 
method used in our numerical simulations, introduced in- 
dependently by West et al. (f8) and by Dommermuth et 
al. (f9| , belongs to this second approach (see also the re- 
cent work byTanaka ^0|). Very recently three new meth- 
ods have been proposed as very promising for simulating 
water waves |2^-|2^] . Results using these new approaches 
on turbulent cascades are still to be completed. 

In this Letter we establish numerically, using a HOS 
method, that nonlinear interactions are sufficient for gen- 
erating power laws in wave spectra; moreover we show 
that the Zakharov-Filoncnko theory is completely con- 
sistent with the primitive equations of motion. We con- 
sider a system of random waves localized in wave num- 
ber space and we show how nonlinearities "adjust" the 
spectrum in agreement with the Zakharov and Filonenko 
prediction. Numerical work in the case of a forced and 
dissipated system has been attempted by Willemsen Q 
using what sometimes are called the "Krasitskii equa- 
tions" (see also Q). In order to avoid the effects of 
external forcing, we considered the case of a freely de- 
caying wave field. If the simulations, as we will see, show 
the formation of a power law then the conjecture that 
this power law is caused by geometrical features related 
to forcing and wave breaking must be excluded, since 
forcing is absent and wave breaking cannot be taken into 
account using the numerical method considered. From 
a physical point of view, the freely decaying case corre- 
sponds to the evolution of a swell wave held. A generic 
wave held is considered at time t = and it is allowed to 
evolve in a natural way using a high order approximation 
of the Euler equations. Since numerical computations are 
limited by the dimension of the grid considered, an arti- 
ficial dissipation is needed at high wave numbers in order 
to prevent accumulation of energy and a break down of 
the numerical code. The fluid is considered inviscid, ir- 
rotational and incompressible. Under these conditions 
the velocity potential ())(x,y,z,t) satisfies the Laplace's 
equation everywhere in the fluid. The boundary condi- 
tions are such that the vertical velocity at the bottom is 
zero and on the free surface the kinematic and dynamic 
boundary conditions are satisfied for the velocity poten- 
tial i/j( x 7 V, t) — 0(^7 2/i il( x i 2/7 *)i t) ( we assume that fluid 



is of infinite depth): 



: \ V ) 2 (1+Vl+V 2 y) 



= (1) 



Vt + tpxVx + ipyVy ~ <f>z | n (l +Vx + Vij) = 0, 



(2) 



The major difficulty for numerical simulations of the sys- 
tem (|f|)-(||) consists in that we have to compute the 
derivatives of <p with respect to z on the surface 77. This 
problem can be overcome if we express the velocity po- 
tential ip(x,y,t) as a Taylor expansion around z = 0. 
Inverting asymptotically the expansion one can express 
4> z \ri as an expansion of derivatives of tp(x,y,t) that can 
then be computed using the Fast Fourier Transform, sim- 
plifying notably the computation. This is nothing other 
than a different way for formulating the HOS method. 
We underline that this is the same approach that has 
originally been used in [Q for deriving analytically the 
equation that is usually known as the "Zakharov equa- 
tion" . The order of the simulation can be decided a pri- 
ori and depends on how many terms are retained in the 
Taylor expansions; in our numerical simulations we con- 
sidered the expansion necessary to take into account four 
wave interactions so that we are consistent with the order 
of the "Zakharov equation" . 

A delicate point in our numerical simulations is related 
to the dissipation of energy at high wave numbers. We re- 
mark that this dissipation is completely artificial since we 
are dealing with a potential flow. Nevertheless we have 
considered the dissipation phenomenon of the wave field 
to be similar to the one that takes place in a turbulent 
flow, i.e. that is mathematically expressed by a Lapla- 
cian that operates on the velocity. As is usually done 
in direct numerical simulations of box turbulent flows, 
in order to increase the inertial range, we have used a 
higher order diffusive term. More explicitly on the right 
hand side of equation (Q)-(||), we have added respectively 
two extra terms: — v(— V 2 ) n ip and — V 2 ) m i], where v 
and [i represent an artificial viscosity coefficient and V 2 
is the horizontal Laplacian. If n and m are greater than 
1 the viscosity is known as "hyperviscosity" . 

It has to be noted that, at first sight, one would use a 
very high power of the Laplacian in order to increase no- 
tably the inertial range, unfortunately very high values of 
to and 71 could bring about the "bottleneck effect" [ p6[ , 
i.e. an accumulation of energy at high wave numbers that 
could distort the power law expected p7J . In our numer- 
ical simulations we used v = /i = 3x 10^ and n = m = 8. 
These values have been selected after some trial and error 
during the development of the numerical code: because 
of our limitated number of grid points, smaller values of 
m and 77, would obscure almost completely the inertial 
range. In our numerical simulations we did not impose 
any a dissipation at low wave numbers. 

In order to prepare the initial wave held it is rea- 
sonable to consider a directional spectrum 5(|k|,0) = 
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P(\k\)G(6). The directional spreading function G{9) 
used here is a cosine-squared function in which only the 
first lobe (relative to the dominant wave direction) is con- 
sidered: 



tions because, as time passes, energy is lost due to vis- 
cosity, thus reducing the significant wave height of the 
wave field and therefore the steepness. 



G(0) 







if -<J <0 <a 
else 



(3) 



(j is a parameter that provides a measure of the direc- 
tional spreading, i.e. as a — * 0, the waves become in- 
creasingly unidirectional. In our numerical simulations 
we selected the value of a — n/2. We tried to avoid the 
complete isotropic case in order to verify if the theory 
still holds for intermediate values of the spreading. At 
the same time the selection of a large value of a was 
motivated by the fact that recently it has been found 
p8| that, for sufficiently narrow angle of spreading, the 
Benjamin-Feir instability can be responsible for the for- 
mation of freak waves. As a consequence the nonlinear 
energy transfer could be slightly altered and some cor- 
rections to the prediction could be necessary (this very 
interesting topic is now under investigation and results 
will be reported in a different paper). P(|k|) is chosen to 
be any localized spectrum. We have performed numerical 
simulations with a gaussian function or with a "chopped 
JONSWAP" spectrum (a JONSWAP spectrum with am- 
plitudes equal to zero for frequencies greater than 1.5 
times the peak frequency) with random phases. For the 
case of the gaussian function, wave numbers lower than 
a selected threshold have been set to zero in order to 
avoid extremely long large waves. The velocity potential 
is then computed from the initial wave field using the 
linear theory. Both gaussian and JONSWAP spectra led 
to the same results in terms of the turbulent cascade. 

Our computation is performed in dimensional units; we 
have selected the initial spectrum centered at 0.1 Hz, i.e. 
we are considering 10 seconds waves. The initial steep- 
ness computed as e = koH s /2 was chosen to be around 

0. 15 (H s was computed as 4 times the standard devi- 
ation of the wave field). The wave field was contained 
in a square grid (the resolution is 256 x 256) of length 
L = 1417.6 meters. The time step considered was 1/50 
the dominant frequency, i.e. At = 0.2 seconds. We have 
performed our numerical simulations on a 400Mh PC. In 
Fig. [I] we show the evolution of the wave power spectrum 
for different real times (t=0, 0.1, 0.5, 1 hours). We see 
that, as expected, the tail of the spectrum starts to grow. 
This process seems to be quite fast: as is shown in the 
figure after a few dominant wave periods some energy is 
already injected into high wave numbers. The process 
of adjusting the power law to the "correct" one becomes 
then very slow, especially for low wave number. This 
could be due to the frozen turbulent phenomenon [p9| , 

1. e. a condition in which the energy fluxes towards high 
wave numbers are reduced because of the discretness of 
the spectrum. Moreover decaying numerical simulations 
are very time consuming with respect to forced simula- 
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FIG. 1. Wave spectra at different times. 



Even though it is not clear from the Log-Log represen- 
tation in Fig. 0, there is a downshifting of the peak of the 
spectrum towards lower wave numbers; as a consequences 
the steepness subsequently decreases over time. The time 
scale of the nonlinear energy transfer becomes larger and 
larger. In Fig we show the power spectrum of the sur- 
face elevation after 4 hours (the steepnes of the wave field 
is e ~ 0.07). In the same plot we show two power laws 
~ k~ 2 - 5 and ~ k~ 3 : the first one seems to better fit the 
data. In order to be completely sure that the numerical 
data are in agreement with the prediction of Zakharov 
and Filonenko, we show in Fig. [| compensated spectra 
with different compensation powers: z = 2.5 seems to 
be the most plausible power law. Thus there seems to 
be ample evidence from our numerical simulations that 
the power law is in sufficiently good agreement with the 
value predicted by the theory. 
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FIG. 2. Wave spectrum at t = 4 hours. A fc~ 2 ' 5 (dot- 
ted-line) and a k~ z (dashed-line) power law are also plotted. 
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FIG. 3. Compensated wave spectra for different values 
of the compensation power: z = 2 (dahsed-line) z = 2.5 
(solid-line) and z = 3 (dotted- line). 

After the pioneering work by Zakharov and Filoncnko 
the kinetic wave theory has developed further, making 
available quantitative predictions for other physical ob- 
servables such as energy fluxes, downshifting of the peak, 
energy dissipation etc. All these quantities will be exam- 
ined and results will be reported in future papers. Other 
questions naturally arise from our results: in HOS sim- 
ulations the order of the computation depends on how 
many terms are retained in the Taylor expansion; do 
higher order terms influence the cascade process? Our 
computation has been performed in a freely decaying 
case; could external forcing (especially if anisotropic) in- 
fluence the power law? And more, what would be the 
influence of the water depth? These are all questions to 
be answered in the near future. 
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